Role of surface waves on the relation between crack speed and the work of fracture 
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We show that the delivery of fracture work to the tip of an advancing planar crack is strongly 
reduced by surface phonon emission, leading to forbidden ranges of crack speed. The emission can 
be interpreted through dispersion of the group velocity, and Rayleigh and Love branches contribute 
as well as other high frequency branches of the surface wave dispersion relations. We also show that 
the energy release rate which enters the Griffith criterion for the crack advance can be described as 
the product of the continuum solution with a function that only depends on the lattice geometry 
and describes the lattice influence on the phonon emission. Simulations are performed using a new 
finite element model for simulating elasticity and fractures. The model, built to allow fast and very 
large three-dimensional simulations, is applied to the simplified case of two dimensional samples. 
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I. INTRODUCTION 

It is well known that fracture propagation is strongly 
governed by properties of the surface created. The clas- 
sical Griffith criterion requires that for a crack to propa- 
gate all the elastic energy supplied to the crack tip must 
at least match the work of surface creation. Crack speed 
is more complicated: one needs to understand the ener- 
getics of mechanisms in competition with surface work, 
and in this paper we focus on phonon emission. 

Phonon emission from the crack tip involves both bulk 
and surface waves. In the continuum limit these waves 
are dispersionless and emission is only expected when 
the crack speed matches any of the speeds vi, Vt and vr 
of longitudinal, transverse and Rayleigh (surface) waves. 
The Rayleigh speed being the lowest, it is often regarded 
as the ultimate theoretical limit for the speed of crack 
propagation. Real materials however have properties 
that do not match the continuum theory. One famous 
example of this is the existence of forbidden crackspeeds, 
which was first explained by M. Marder et alBuB taking 
into account the discrete nature of matter. Their qual- 
itative explanation of the velocity gap where no crack 
can propagate in a periodic lattice at speeds lower than 
roughly 1/3 of the transverse waves speed, was related 
to lattice oscillations breaking bonds before the expected 
time compatible with the crack speed. This explanation 
involves the use of a "most stretched bond" breakage rule. 
Here we propose a more cautious but general view using 
an approach based on the phonon band structure and 
the energy release rate. The phonon band structure can 
lead to resonant emissions which influence the fraction 
of energy radiated from the crack tip. The existence of 
forbidden crack speeds then follows. 

When dealing with dispersion relations the reference is 
usually to bulk waves or, in the case of surface waves, to 
the Rayleigh branch responsible for the Rayleigh speed. 
Geophysicists and researchers in earthquakes are also fa- 
miliaE-jwith another kind of surface waves known as Love 
wavesatl, polarized in the direction normal to that of 
propagation and parallel to the crack plane. Love waves 



only exist in the continuum limit when there is a gra- 
dation of elastic properties near the surface, pertinent 
to seismology but not usually considered in fracture me- 
chanics. However, due to the discrete nature of matter at 
the atomic level dispersive media can also support Love 
and other wave branches so that they can enter fracture 
problems. It is natural to include the whole complex- 
ity of the surface dispersion relations in the description 
of the material to understand its role, and how and if 
this complexity can change the theoretical description of 
fracture dynamics. 

The form of dispersion relations is the fingerprint of the 
discrete spatial organization of matter, beyond the con- 
tinuum limit. In any kind of simulation it is necessary 
to deal with discrxlization in space. In molecular dy- 
namic simulationsmffi] this corresponds to an atom istic 
description; in mesoscopic.mjacLels as lattice-likeEJii3l!3'L3 
or finite element modelsc3li3'tZl, discretization has to be 
imposed in space in order to solve differential equations 
of the dynamics. This has an unavoidable effect on the 
dispersion relations that necessarily reflect the underly- 
ing lattice. 

In this paper we present a new model developed to sim- 
ulate an elastic continuum material using finite tetrahe- 
dral elements. Space is discretized on an fee grid, and this 
is reflected in the dispersion relations of the simulated 
material. This feature could be considered a drawback, 
but can be exploited to model the influence of similar 
dispersions in real fee structure materials. The simplic- 
ity of the model allows us an extensive understanding of 
both simulations and theory. The model, built to allow 
fast and very large three dimensional simulations, is used 
here in the simplified case of two dimensional samples, 
and phonon emission from the crack surface of planar 
cracks advancing at fixed speeds is analysed. An analy- 
sis of the influence of this emission on the crack dynamics 
is the aim of this work. 

New and central to our approach is that we fix the 
crack properties (speed and shape) and measure some of 
the mechanical properties (stress, energy release rate). 
This can be done with success in two-dimensions and 
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gives new insight into the behaviour of some of these 
quantities. In particular, we show that the energy re- 
lease rate is strongly influenced by the surface phonon 
emission at the crack tip through resonant emission at 
particular crack speeds, leading to the existence of nar- 
row bands of permitted crack speed. We also show that 
the energy release rate which enters the Griffith criterion 
can be expressed as the continuum solution, multiplied 
by a microscopic function that only depends on the lat- 
tice geometry and that describes the lattice influence on 
the phonon emission. 

The paper is organized as follows: in section II a full 
description of the model is provided, including the de- 
scription of the bulk dispersion relations. In section III, 
a preliminary analysis of simulations of planar cracks at 
constant speed is given, and in section IV we relate the 
phonon emission to the dispersion relations of surface 
waves and we show the phenomenon of resonances. In 
section V we analyse how the energy release rate is influ- 
enced by the phonon emission and we derive the existence 
of bands of permitted crack speeds. Finally, conclusions 
and comparison with other works are given in section VI. 



II. SIMULATION MODEL 

The model is an application of the finite element 
method with the aim of simulating linear elasticity and 
fractures. In continuum linear elasticity!!! the Lagrangian 
has the form: 



L = 



pu_ 

2 



"2 VU 



dV 



where the stress tensor is given b; 



,yi: 



a= ATr(Vu) 1 + /i[Vu + (Vu) T ]. 



(2.1) 



(2.2) 



The unsymmetrized strain tensor Vu from which the 
stress tensor, Lagrangian and hence equations of motion 
follow is central to our finite clement scheme below. 

The key driving features of the model are that the elas- 
tic response is as local as possible so that non-linearity 
and particularly rupture can be incorporated, whilst ac- 
cidental soft modes are avoided and all mode frequencies 
are strictly upper bounded for stability of timestepping. 




FIG. 1: Unit cell for the fee lattice used in our simulations. 
This contains eight tetrahedra, each connecting four lattice 
sites. Displacement and momentum are defined on the sites, 
whilst strain and stress are defined on the tetrahedra. 



where the v index spans over all lattice points and the t 
index spans all the tetrahedral elements, fl' is the volume 
of the system per tetrahedron, 1 /8 of the volume of the 
fee conventional unit cell shown in fig. |l|. 

The stress tensor is eva luated at the center of each 
tetrahedron through (2.2), given the unsymmetrized 
strain tensor at the center of the t-th tetrahedron is: 



(Vu) 



o^y^ r tv u v 



(2.4) 



with the condition a J2 V r tv r tv = 1, where r tv is the vec- 
tor joining the center of the t-th tetrahedron with neigh- 
bouring vertex site v, and u„ is the displacement field 
on site v. Finally, the equation of motion takes then the 
form: 



v = f v = -aft 1 ^2 r tv 't 



raw.. 



where a t is the stress tensor of the tetrahedron. 



B. Dispersion relations for bulk waves 



(2.5) 



A. Definition of the elastodynamic model 

To evaluate the unsymmetrized strain tensor, space 
is discretized using an fee lattice and lattice points are 
connected using tetrahedral elements. Each element con- 
nects four points (see fig. ^), which is the minimum for 
a full local gradient calculation in three dimensions. On 
each lattice point, the displacement field u is defined, and 
the Lagrangian in its discretized form is given by: 



L 



(2.3) 



The dispersion relations of our model are relatively 
straightforward to calculate. These confirm that in the 
long wavelength limit the continuum case is retrieved, 
and particularly that the use of tetrahedral elements 
avoids the presence of soft modes at the border of the 
Brillouin zone, a problem arising when using simple cu- 
bic schemes. 

Introducing a spatial Fourier transform in equation 
~Th, the gradient operator can be written as: 



FT[V] = a^r„e' 



kr. 



(2.6) 
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FIG. 2: Acoustic dispersion curves for the tetrahedron model along different directions in reciprocal space: there are no soft 
modes at the border of the Brillouin Zone. These plots are shown for Poisson's ratio v = 1/3 and using units where a = 1 for 
the side of the tetrahedral cell and v t = y2 for the speed of transverse acoustic waves. 



There are two kinds of tetrahedra (see fig. |j) which 
can be readily seen to be mirror images of each other, 
and their corresponding gradient operators are related, 
in the Fo urie r transform, by Vi = — V|. In terms of 
these eq. (2.5) can be expressed as: 



iu = n' 



V 2 - O"! +V r CT; 



where 



<Ti= A / Vi • u + /i 



Viu + uV] 



(2.7) 



(2.8) 



and £2 is obtained by changing the index from 1 to 2 in 
eq. (P). 



Equation (2/7) is therefore given, after Fourier trans- 
forming, by 

mw 2 u k = fi' j/i V • V* + V* • V u k + 

+ (A + ^)[w* + V*vl u k l (2.9) 



where V = Vj = — V2. This equation can be solved ex- 
actly for eigenvalues (see appendix^]), giving the results 



shown graphically in figure g. The dispersion curves are 
well behaved, and reflect the underlying fee structure: 
even if we are trying to simulate the continuum limit, 
dispersion relations are the image of the spatial struc- 
ture of the underlying lattice. This cannot be avoided in 
our model just as in real materials. 

The dispersion relations found vindicate the use of 
tetrahedral elements. Using a simpler cubic lattice and 
cubic finite elements leads to unphysical soft modes, 
where the frequency goes to zero whenever two of the 
three components of the k vector approach the border of 
the Brillouin zone. The use of tetrahedra removes this 
problem since there is no way to deform a tetrahedron 
giving a null contribution to the strain tensor, whilst al- 
ternative strategies involving less local spatial derivatives 
make fracture properties harder to implement and con- 
trol. 

The value of the maximum frequency is set by the lon- 
gitudinal branch at th e Brillouin Zone boundary. We 



will see in section II D that such a value remains an up- 



per bound even when considering the finiteness of the 
sample or the presence of fractures. This value provides 
the minimum period of site motion. 

Due to the significant role of surface waves in this pa- 
per, a study on their dispersion relations will be described 
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later in section IV 



C. Discretizing time 

Time can be discretized by introducing a finite 
timestep h. The integration scheme used in our simula- 
tions is the Leap-Frog schemet£l, in which displacements 
u„ and momenta p„ are evaluated alternately at subse- 
quent steps: 



u v (t) 



= u v (t- 2h) + 2h 



Pv{t - h) 



p v {t + h) = p v (t-h)+2M v (t). 

Fourier transforming these equations and comparing 
with the limit h — > leads to the dispersion relation: 



sin (uh) = h uj^(k) 



(2.10) 



where &t>o(k) is the value of the frequency of mode k in 
the continuous time [h — > 0) limit. From this equation 
it follows that we require /iwo(k) < 1 for mode k to be 
stable in the discrete timestep simulations, and the global 
stability limit is 



h< 



max wo(k) 

k 



It is interesting to apply this to the simplified case of 
one dimension. In this case the continuous time disper- 
sion relations have the simple form 



co 2 (k) 



sin 



(f) 



(2.11) 



with: u>max = (2/a)\/ (A + 2p)j p. Comparing ( 2.10 ) and 
( 2.11 ), it is easy to observe that the stability limit for 
the timestep is h — l/w m ax- Moreover, if we use exactly 
this value for the timestep, we obtain the linear disper- 
sion relation u) = ck with c = aLu max /2. Thus using the 
maximum allowed timestep, the relation between the fre- 
quency and the wave vector is the same as in continuum 
elasticity. This notable property is valid only in one di- 
mension and is lost in two and three dimensions where 
the dispersion is more rich. Nevertheless, this property 
shows that having discretized space, it not necessarily an 
optimal approximation to the continuum to use a very 
small timestep as in standard Molecular Dynamics sim- 
ulations. 

In three dimensions w m ax is given by the maximum 
frequency over the entire first Brillouin zone. Using the 
maximum allowed timestep dispersion relations are closer 
but not equivalent to the continuum limit. 

The Leap-Frog scheme used here to discretize time is 
relatively robust with respect to energy conservation due 
to its time reveribility and symplectic properties. E3 There 
is an exactly conserved-fe-dependent hamiltonian Hh very 
close to the naive oneEa, with for example: 



H h 



1 



K+ + K\ + U + 0(h 2 ) 



(2.12) 



where K + and K_ are the kinetic energies advanced and 
retarded by ±h relative to the timestep where the poten- 
tial energy U is evaluated. 

In these simulations the value chosen for the timestep 
has been h = 0.1 /imax- Using such timestep, dispersion 
relation are close to the continnum time dispersion rela- 
tions. Some trial simulations with h = 0.01 /i max have 
been performed showing no substantial difference with 
the results here reported. 



D. Nonlinear elasticity and breakage 

The model up to now has been described as a model 
for linear elasticity. However, since the link between 
the stress and t he s train tensor is given exclusively 
through equation (2.2), we can easily generalize the equa- 
tion to include any kind of elastic response as well as 
anisotropies. For instance: 



cr=f (Vu+ (Vu) rj 



(2.13) 



where each tensor component fij represents any kind of 
nonlinear function. In this paper we study only the sim- 
pler scenario of linear elasticity augmented by breakage 
of tetrahedra corresponding to the advance of the crack. 
We brake tetrahedra by abruptly setting its elastic con- 
stants A and p, to zero and it is the resulting recoil of the 
neighbouring sites which excites phonon emission. 

Breaking tetrahedra does not compromise the maxi- 
mum allowed timestep, because we can readily show that 
the maximum vibrati ona l frequency cannot increase. To 
see this consider eq. (2.5) rewritten in the form 



pu v = 



dU 

8Uy ' 



From (^), we clearly have ^2 V u v ■ dU /du v = 2U, where 
hi is the total potential energy. Now consider a normal 
mode obeying: 



dU 
du v 



from which it follows that 



put 



2U 
J2v l u " 



Breaking a tetrahedron removes strictly non-negative 
terms from the numerator of the right hand side, and 
this can only lower the maximum of to. 



III. SIMULATIONS OF PLANAR CRACKS AT 
CONSTANT SPEED 

An advancing crack is simulated through the break- 
ing of tetrahedra. To decide when to break we need a 
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FIG. 3: Set-up for simulations. A fixed displacement is ap- 
plied on the top and bottom boundaries, and the crack ad- 
vances at constant speed up to the middle of the sample. 



breakage rule. One possibility is to break tetrahedra as 
soon as they overcome some fracture criterion, for exam- 
ple a critical value of principal stress or elastic energy. In 
this work we instead follow a complementary approach, 
advancing a straight crack at constant speed and measur- 
ing the fracture criteria achieved such as work of fracture. 
This provides new insight into the relationship between 
energy release rate and crack speed, including any inter- 
vals of forbidden velocities. 

For computational convenience we focussed on crack 
speeds commensurate with our timestep, that is 



1 

2hn 



for n = 1,2,3,... . Any desired crack speed can be ob- 
tained from this sequence by modest adjustment of the 
timestep h. 



A. Stress field of advancing cracks 

The model presented above is fully three dimensional. 
However for a strictly planar type-I crack, commensurate 
with our lattice, it is readily shown that the resulting dy- 
namical solution is strictly planar with the displacement 
all in that plane. We have exploited this and shrunk our 
system to one fee cell deep in the third dimension, with 
periodic boundary conditions. 

In our simulations the Lame coefficients are set so that 
the Poisson ratio is v = 1/3. The boundary conditions 
correspond to imposing fixed normal counter displace- 
ment at the top and bottom boundaries, i. e. u x = 0, 



±o at y — ±i/2, and a xx = a x 



on the right 



and left faces (see fig. |3|). The initial condition is the 
static elastic solution, found by relaxing the lattice to 
its configuration of minimum energy by adding dissipa- 
tion as discussed in appendix [b|. Because of the linearity 
of the equations, the magnitude of the imposed bound- 
ary displacement does not influence the results: for the 
following results the starting displacement is 1 % of the 
sample height. The advancing fracture is simulated until 



it reaches the middle of the sample, then a snapshot of 
the stress field is taken. 

In figure ^. (a) the trace of the stress field for a typical 
simulation is shown: strong emission of surface waves is 
visible on the crack surface, as also shown in figure ^.(b). 
In particular, from this last figure it can be seen that the 
amplitude of these oscillations is smaller but of the same 
order of magnitude as the stress level at the crack tip. 



B. Description of the basic phenomenon 

When a crack is advancing new surface is being cre- 
ated, and waves are emitted both into the bulk and along 
the crack surface. In a system of reference comoving with 
the crack tip, the frequency u)' of emitted waves has to 
match the temporal frequencies with which lattice struc- 
ture presents itself to the crack tip, leading to the selec- 
tion rule: 

w = w — k ■ v = g ■ v 

where g is any reciprocal lattice vector. This is equivalent 
to a simple matching of phase velocity: 

w = kv (3.1) 

where the wave vector k is viewed in the extended zone 
scheme as sketched in fig. |[ 

From a graphical point of view, a crack advancing at 
a given speed can be represented as a straight line in 
the {k x , uS) plane, its slope given by the crack speed. The 
crack emits waves at the frequencies and wavelengths cor- 
responding to the intercepts of this line with the disper- 
sion relations (see fig. ||) in the exte nde d zone scheme. 

For some crack speeds condition (3.1) may be fulfilled 
not just for one k point, but over a neighbourhood of k. 
This can readily be shown to reduce to: 

dui 

v G = ^ = v. (3.2) 

This matching of phase velocity and group velocity 
(sketched as well in fig. |^) is the condition of resonant 
emission. We can expect for these crack speeds a sharp 
increase in the intensity of the phonon emission. 

To take this discussion further we need to describe the 
dispersion relations of surface waves for our model, as 
detailed in the following section. 



IV. DISPERSION RELATIONS FOR SURFACE 
WAVES 

The description of the modes in which a lattice with 
an underlying fee symmetry can vibrate is more complex 
than the one given by continuum elasticity. Frequency as 
a function of wave vector is no longer linear, and becomes 
periodic in the extended zone scheme. Moreover, new 
branches of surface wave appear besides the continuum 
Rayleigh branch. This richer band structure is crucial in 
understanding features on resonant phonon emission. 
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FIG. 4: (a) Absolute value of the trace of the stress field for a crack advancing at half of the transverse wave speed, with 
Poisson ratio v = 1/3. This is the central 200 x 200 region from a sample of 400 x 400 tetrahedra. Intensity spans from white 
(zero) to black (highest values), (b) Trace of the stress field measured along the fracture plane, through the crack tip whose 
singularity is clearly visible, and onwards ahead of the crack tip. 



A. Theoretical surface dispersion relations 

Due to the simplicity of the model it is not difficult 
to evaluate the dispersion relations of its surface waves. 
These are solutions of w 2 Uk + V- o— with boundary 
condition o -n = 0, where n is a unit vector normal 
to the surface. Solutions are linear combination of bulk 
modes characterized by the same frequency ui and surface 
wave vector components, and different complex k- n that 
give decay towards the interior of the sample. 

Calculations of the surface modes are reported in ap- 
pendix ^|. Results are shown in fig. |g.(a). 

The figure shows quite a rich diagram which includes 
the Love branch and some high-frequency branches. It is 
clear that all the branches could take part in the process 
of phonon emission. 



B. Direct measurement of surface waves 

The dispersion relations can be directly measured by 
exciting the top free surface of a simulated elastic sam- 
ple. Sites of the top face are displaced according to one 
specified wave vector k x . The surface is let free to os- 
cillate and the frequencies of the corresponding excited 
modes are measured. Changing the k x vector at k z = 0, 
the entire first Brillouin zone can be spanned and a di- 
rect measurement of the dispersion of surface waves can 
be obtained. Results are shown in figure |[(b). 

It is evident that all the surface excitation is associated 
with the theoretical branches. No emission correspond- 
ing to the Love branch is visible due to its polarization 
perpendicular to our two dimensional plane. We have not 



identified the reason why no emission is observed corre- 
sponding to the Riv branch. 

Our measurements do show continuation of the theo- 
retical branches beyond their end points. We have veri- 
fied from calculations of appendix |C] that these continu- 
ations correspond to frequencies with a small imaginary 
part. 
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FIG. 5: Schematic description of the basic matching phe- 
nomena. A crack advancing at some speed v matches wave 
phase velocity at points A, B, C and D, so that correspond- 
ing phonons can be emitted from the crack tip. The figure 
shows the special case in which there is a resonant match at 
D: the crack speed also matches the wave's group velocity 
which should lead to strong emission. 
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FIG. 6: (a) Points represent theoretical dispersion relations for surface waves in the first Brillouin zone, for our model with a 
(0, 1) surface normal and v = 1/3, with units as in figure^. Lines show the dispersion relations for bulk modes. All surface modes 
have Im(fcy) > to be damped for y — > oo, and k z — corresponding to the two dimensional constraint in our simulations. 
The Rayleigh and Love branches are labelled as Rj and L. The other high frequency branches are not interpretable in the 
continuum limit, (b) Dispersion relations measured by excitation of a free surface, compared with the theoretical branches. 
More than one frequency can be excited for a given k x . The size of circles reflects the intensity of the signals. Only three of 
the five branches are excited: the Love branch is not visible due to its polarization, however it is not clear the reason of the 
lack of evidence of the Rjy branch. 



C. Emission from advancing cracks 

Let us consider a snapshot at time t of the crack sur- 
face. Depending on the crack speed, some of the vibra- 
tional modes will be excited. Through a Fourier trans- 
form of the surface profile th e ex cited fc-vector can be 
obtained. Assuming that eq. (3.1) applies we can com- 
pare observed values of k x at crack speed v with k x ver- 
sus Lu/k x from the theoretical dispersion relations. The 
one delicate point is that k x is indeterminate up to mul- 
tiples of a reciprocal lattice vector g. Figure shows 
the agreement between measured emission and the the- 
oretical curves, with assignment of multiples of g to the 
measured k x as the only element of fitting. 



D. Resonances 



Resonances appear as soon as eq.(|3.2|) is fulfilled: this 
matching of phase velocity and group velocity can lead 
to resonant emission as sketched in fig. |[ Plotting the 
phase velocity co/k x versus k x , the resonant condition is 
equivalent to a plateau in uj/k x . Predicted resonances 
can therefore be easily read from our computed disper- 
sion relations as shown in fig. ^ Due to the periodic 
band structure there is an infinite number of possible res- 
onances corresponding to the infinite number of Brillouin 
zones. However we observe (not surprisingly) that signal 
corresponding to resonances in higher Brillouin zones is 
weak, so our discussion in this paper focuses on the first 
two Brillouin zones only. 

The presence of resonances can be verified by analysing 
the total intensity of emitted waves as a function of crack 



speed. Figure g shows the total wave intensity integrated 
over just the crack surface, whilst the discussion of the 
full radiated wave power will be given in section 

Corresponding to the speeds where resonances are ex- 
pected there is an effective peak in the wave intensity. 
These peaks are slightly shifted towards lower speeds, 
which can be interpreted in terms of the concavity of 
dispersion relations at the resonance. When the crack 
advances at a speed slightly below the resonance, more 
modes can be excited, giving a maximum in intensity 
slightly shifted. The abrupt drop of intensity at higher 
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FIG. 7: In this plot crack speeds correspond to horizontal 
lines. Black dots correspond to the emission measured from 
the crack surface. For simplicity only emission from the first 
and second Brillouin zone (folded back into the first zone) 
has been reported. All dots lay on the theoretical dispersion 
relations or on their extensions in the second Brillouin zone. 
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crack speeds (approaching the Rayleigh speed) appears 
to result from the fact that the penetration length of such 
waves is of order their wavelength, and this in turn be- 
comes of the order of the sample size as the Rayleigh 
speed is approached. At high speeds, the total intensity 
of surface waves cannot be retrieved by analysing just the 
crack profile. 



V. ENERGY RELEASE RATE 

The presence of resonances suggests that an effect on 
the crack dynamics has to be expected. We will now show 
how the most commonly and relevant measured quantity 
of crack dynamics, namely the energy release rate, is af- 
fected by this phenomenon. 

In the simulations here reported, a planar crack is left 
to advance at constant speed from a short starting notch. 
At early times the crack is almost equivalent to a crack 
advancing in an unbounded medium, since the sample 
boundaries are distant from the crack tip. When the 
crack becomes large compared to the linear dimensions 
of the sample the crack starts to "feel" the presence of 
the boundaries and the correct description is that of a 
crack advancing in a strip. In the continuum and steady 
limit the short crack regime is characterized by a linear 



X(v) 



8 


5 




! P 


) a' a 


. 


2 




0.4 0.6 


0.8 



FIG. 9: The total wave intensity integrated over the fracture 
surface, as a function of crack speed. Peaks in the emission 
correspond to the resonances shown in figure ^| as indicated 
by vertical lines. A systematic shift towards the lower crack 
speeds is visible and is related to the convexity of the disper- 
sion relations (see text). These results were obtained from 
simulations of samples of 800 x 800 tetrahedra. 

increase of the energy release rate G with time, whilst in 
the long crack regime the G function is time-independent. 
Both regimes are sensitive to the crack speed. 
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FIG. 8: The theoretical phase velocity uj/kx versus k x , with 
horizontal lines showing where resonant phonon emission can 
be expected. The lines show resonances due to the first and 
the second Brillouin zone only. Dashed lines correspond to 
resonances due to bulk dispersion relations. 



1. Energy release rate in a discrete sample 

The discrete nature of matter is reflected in the de- 
pendence of the energy release rate G on crack speed. 
We can describe the macroscopic energy release rate per 
unit distance of crack advance Gm as the sum of two 
contributions: 



G M (v,t) = G hT (v,t) + G ph {v,t). 



(5.1) 



Gm{v, t) is the solution of the continuum limit which gov- 
erns the macroscopic delivery of energy towards the crack 
tip, for which we have theoretical expressions available. 
G^v, t) and Gp n (t>, t) are respectively the breakage en- 
ergy release rate and the phonon energy release rate. Our 
strategy below is to directly measure the breakage energy 
release rate from the potential energy lost when tetrahe- 
dra are broken. 

The macroscopic energy release rate Gm(v, t) is largely 
determined by the macroscopic conditions and the length 
of the crack. In the long crack limit in the case of a strip 
of height 21 and fixed displacement 8 at each boundary, 
the macroscopic energy release rate corresponds to the 
am^tuit of elastic energy stored far ahead of the crack 



ti 



tip has 



From eq.(2.2) the stress field ahead of the crack 



(A + 2fi)e 



(5.2) 



where e = 8/1 is the imposed strain. Hence the time 
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FIG. 10: The measured energy delivered into bond breaking 
per unit length of crack, as a function of time for different 
crack speeds. From these data we measured the early time 
slopes G'^ r (v, 0) corresponding to the short crack regime and 
(from the slope of the time integrated plots inset) the long 
crack plateau values G , j 3r («). The (non-trivial) sequence of 
the curves is more readily appreciated from fig. pj[ 



independent macroscopic energy release rate is: 

Gm - (A + 2fi)e 2 l (5.3) 

independent of the crack speed v. 

The translation of the Griffith criterion in the discrete 
case is that the crack will advance as soon as the energy 
stored in the tetrahedron to break, G^j., is greater than 
a threshold value connected with the toughness of the 
material. The behaviour of G^j. with the crack speed is 
therefore crucial for understanding the crack dynamics. 

The breakage energy release rate clearly depends on 
the external macroscopic conditions as well as the dis- 
creteness of the model. However it can be expressed as: 



G hT (v,t)=E(v)G M (v,t) 



(5.4) 



where we have introduced the efficiency E(v). In the 
long crack limit for the fixed grip set-up used in our sim- 
ulations this gives Gv jr (w,t) = E(v)Gfj where the effi- 
ciency is the sole source of velocity dependence. The ef- 
fect of local discreteness thus separates from the effect of 
macroscopic external conditions: all the dependence on 
the crack speed is hidden in the efficiency function E(v) 
which is local to the crack tip region and independent 
of the macroscopic regime as will be shown below. The 
meaning of the efficiency E(v) is as follows: when E(v) 
is close to zero, the energy delivered to the crack tip is 
mostly spent in phonon emission so that the mechanism 
is not sufficient for the crack to advance. When the E{v) 
is close to one, all the energy delivered is used to break 
tetrahedra, and the crack can advance promptly. We will 
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FIG. 11: Dependence of breakage energy release rate on crack 
speed in the long crack regime. Crack propagation with a 
steady speed should only be stable where this function is de- 
creasing. Vertical lines show the theoretical position of reso- 
nant surface wave emission: the two dashed lines correspond 
to resonant bulk wave emission. The figure equivalently shows 
the behaviour of the efficiency E(v) with the crack speed, as 
for long cracks and fixed grip conditions the two are strictly 
proportional. 



see that the dependence on v is crucial in the determina- 
tion of bands of permitted crack speeds, not described in 
the continuum elastic theory. The correct derivation of 
these bands will be given below. 



2. Measuring the energy release rate 

The breakage energy release rate G] 3r (u, t) in our simu- 
lations corresponds, for a given crack moving at speed v, 
to the value of the elastic energy that disappears from the 
system with each broken tetrahedron at time t. Figure 
[To| shows measurements from a set of simulations involv- 
ing samples of of height 120 tetrahedra and up to 3500 
tetrahedra long for increasing crack speeds. We observe 
G]3 r (ti, t) to grow linearly with time initially: correspond- 
ing measured slopes GL.(z;,0) are discussed below. The 
long crack regime exhibits fluctuations due to waves re- 
flected from the sample boundaries, but the average mag- 
nitude Gbj.^) of the breakage energy release rate can be 
clearly retrieved. 

The Gj3 r (u) in the long crack limit is shown as a func- 
tion of crack speed in figure [ll]. Vertical lines indicate 
the speeds at which resonances are expected to be seen 
according to figure |] The two lines labelled as a' and f3' 
correspond to resonances due to the dispersion relations 
of bulk waves. Corresponding to each line there is a clear 
decrease of the energy release rate, meaning that more 
energy is emitted as radiation. 

The crucial feature of figure [ll] is that it reveals how 
the energy available for breaking bonds responds to the 
speed of the crack, for a given macroscopic energy release 
rate. This same figure can be read in reverse: given some 
threshold value for G^ T (v) corresponding to a given frac- 
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ture toughness of the material, the possible crack speeds 
are obtained from the graph. 

It can be further argued on grounds of stability that 
only the speed ranges where G^ r (v) decreases with v 
are allowed as steady crack speeds. In the counter case 
where G\y T (v) increases with v, a prospective overshoot 
in v leads to excessive bond breakage energy and hence 
acceleration, and an undershoot to insufficient breakage 
energy and the crack must slow further. Thus stable 
steady crack propagation is confined to narrow intervals 
associated with resonances and the high speed regime. 

It is worth pointing out that most of the features shown 
in figure [ll] are due to high frequency branches either in 
the first or second Brillouin Zone. The first resonance 
due to the Rayleigh branch (apart from the resonance 
at the Rayleigh speed) is the marginally visible e reso- 
nance. This shows how important it is to include the full 
complexity of the band structure within the analysis. 

The data in fig. [ll] are at constant macroscopic en- 
ergy release rate Gjg, so due to (5.4) they also show the 
behaviour of the efficiency E(v) with the crack speed. 
However, if the efficiency E(v) is governed only by lo- 
cal phenomena, then its behaviour with the crack speed 
should be independent of the crack propagating in steady 
state or transient regime, provided the crack speed v is 
fixed. From the knowledge of E(v) then we could find out 
the dependence of G-^ v {v,t) on the crack speed for any 
macroscopic set up, provided we know the macroscopic 
continuum solution. 

To verify the independence of E(v) on the dynamic 
regime, we also analyzed the short crack limit. In this 
case, the crack can be seen as advancing in an unbounded 
medium as the sample boundaries are distant from the 
crack tip. For a planar type-I crack propagating at steady 
speed v in an unbounded continuum medium the macro- 
scopic transfer of energy to the crack tip (per unit dis- 
tance advanced) is given by: 



where R(q) = {b 2 - 2q 2 ) 2 + <iq 2 ^fa? - q 2 yjb 2 - q 2 and 



G%{v,t) 



iron 
2c 2 iiD 



K 2 (v,t). 



(5.5) 



Here a* , 



1 



v 2 /v 2 l and D 



4a t ai~{l + a 2 ) 2 . The 



constants vi and v t represent the longitudinal and trans- 
verse sound speeds, and Kj{v,t) is the stress intensity 
factor at time t for the given crack speed. 

A functional fopn_for Kj(v,t) is available for the 
Broberg problemotBO of a crack expanding from zero 
initial length in a uniform tension field and infinite 
medium. In particular it can be written in the form: 

Ki(v,t) = S(w)croo\/7rwi 

where Coo is the traction applied on the crack faces. The 
energy release rate increases linearly with time. As figure 
|l0| shows, we obtain the same behaviour for G^ r (ii, t) in 
our simulations. In the case of symmetric growth T,(v) 
has the form: 



£(«) = - 



I{b/q)R{q) 



R(i v ) 



(q 2 + t] 2 ) 3 / 2 a? + if 



dr) 



with a 



q 



We thus have a closed 



form for the early time slope G' M (v,0). 

Due to these results, to prove that E(v) is independent 
of the dynamic conditions we just have to check that 



E(v,0) 



G' hT (v,0) 



G' M (v,0) 



(5.6) 



b 2 q^J q 2 — a 2 



matches E(v) already measured from the long crack limit. 
Simulations of symmetric crack growth can be performed 
simply by fixing the longitudinal displacement at the left 
and right boundaries: though this would correspond to 
having a periodic system of symmetric growing cracks, 
the interaction between cracks is weak due to the strip 
geometry. G'^ (v,0) then corresponds to the measured 
slopes of the early regime, and G' M (v, 0) is the time 
derivative of the macroscopic energy release rate ( p.5| ). 

Results of the comparison are reported in table |TBe- 
cause the product vt corresponds to the crack length, the 
case v = was reconstructed by analysing several sim- 
ulations of static cracks with different crack lengths and 
retrieving the resulting overall slope. For low and high 
crack speeds, problems arise when measuring the short 
crack regime: for low crack speeds data is noisy and there 
are few points to average; for high crack speeds the long 
crack limit is achieved soon after the crack starts mov- 
ing, so that the short crack regime is properly defined 
only when the sample height is very large with respect to 
the size of a tetrahedron. This is still not the case even 
when dealing with our 120 tetrahedra wide sample. Bet- 
ter results for the short crack regime were obtained from 
simulations of 800 x 800 tetrahedra as the third column 
of the table shows. 

These results show that the efficiency E(v) is inde- 
pendent of the particular dynamic regime. The abso- 
lute values are also of interest. The v = limit shows 
that when the tetrahedra at the tip of a static crack are 
broken, about 75 % of the strain energy released comes 
from relaxation in other tetrahedra which is radiated in 
waves. The maximum efficiency (about 60 %) occurs at 
v ~ 0.6 Vt- This happens to be just below the Yoffe speed 
~ 0.63 ft — 2/3wfl, but as our measurements exclude the 
possibility of crack branching we presume this to be a 
coincidence. 



VI. CONCLUSIONS 

We have presented a new finite element model for lin- 
ear elastic fracture mechanics, which has proved surpris- 
ingly and revealingly rich in its behaviour even in the 
two dimensional case. The model was designed to enable 
fast numerical simulations of large systems particularly 



11 





21 = 120 


21 = 800 


21 = 120 




v/v t 


G' M (v,0) 




E(v) G5(«) 


E(v,0)/E(v) 


0.0 


0.257 ± 0.006 




0.2591 ± 0.0008 




0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 


0.19 ±0.04 
0.52 ±0.03 
0.40 ±0.02 
0.44 ±0.02 
0.48 ±0.02 
0.59 ±0.02 
0.55 ±0.04 
0.2 ±0.1 


0.22 ±0.04 
0.52 ±0.01 
0.390 ± 0.006 
0.432 ± 0.004 
0.462 ± 0.004 
0.576 ± 0.006 
0.514 ±0.006 

0.33 ±0.02 
0.079 ± 0.007 


0.2408 ± 0.0004 
0.528 ± 0.003 
0.379 ± 0.001 
0.4315 ± 0.0002 
0.4636 ± 0.0001 
0.5767 ± 0.0004 
0.5201 ± 0.0001 
0.3202 ± 0.0001 
0.0890 ± 0.0005 


0.9 ±0.2 
0.98 ± 0.02 
1.03 ±0.02 
1.00 ±0.01 
0.997 ±0.009 
1.00 ±0.01 
0.99 ±0.01 
1.03 ±0.06 
0.89 ± 0.08 



TABLE I: The efficiency of energy delivery into bond breaking, compared between the different dynamical regimes of short and 
long cracks, as a function of crack speed. The second and third columns show the efficiency computed from early time slopes of 
the energy into bond breakage compared to that macroscopically delivered, for two different sample heights. Results for v = 
were obtained from a set of simulations of static fracture with different notch lengths. Results for low and high crack speeds 
are influenced by problems in the measurements in the short crack limit: this effect is reduced increasing the size of the sample 
as shown in the third column. The final column shows the ratio of these short crack regime efficiencies to those computed from 
the long crack regime shown in the fourth column: these results clearly confirm that the efficiency is insensitive to dynamical 
regime, consistent with our assertion that it is a property local to the crack tip. 



with three dimensions in mind, but the main strengths 
which we have exploited in this paper are direct control 
of the local physics which it offers, combined with bulk 
and surface dispersion relations amenable to simple theo- 
retical computation. Our results suggest that inclusion of 
known phonon dispersion can be crucial to understanding 
the speed of fracture propagation. 

Our results relate to ideally brittle materials, in that 
we have included no significant mechanism of local dis- 
sipation and most particularly no plastic deformation 
mechanism. Linear damping is readily included, and in- 
deed could be exploited to mitigate the effects of waves 
reflected back from the sample edges. Spontaneous crack 
roughness and branching will be addressed in a following 
paper. 

The crucial mechanism which our results incorporate 
beyond continuum fracture mechanics is the radiation of 
phonons from the crack tip. This we show leads to a sig- 
nificant speed dependence in the fraction of macroscopic 
strain energy available as work to create new surface. 
For a static crack in our model, this efficiency is only 
25 %, thus modifying the most naive Griffith criterion 
for crack propagation by a factor of 4. The general rise 
in bond breaking efficiency with speed towards a global 
maximum for typical fast crack speeds rules out steady 
crack propagation at most lower speeds, and remains to 
be understood more quantitatively. 

We have been able to interpret fine structure in the 
bond breaking efficiency associated with resonant condi- 
tions for surface wave emission. This leads to islands of 



stable crack velocity, which can only arise at lower speeds 
due to dispersion and hence due to the discrete properties 
of matter. 

A related explanation of the presence of sets of forbid- 
denx£|ack speeds has already been given by M. Marder et 
aim by using a mechanistic description of bond break- 
age in terms of most stretched ones. According to this 
approach cracks cannot advance at speeds below some 
threshold, since due to lattice oscillations bonds would 
break before the expected time compatible with the crack 
speed. A velocity gap then appears. The two approaches 
seem to be very close indeed, since simulations from Mur- 
der et al. are based on periodic lattices which should 
show characteristic dispersion relations. The oscillations 
between neighbouring sites which stretch the bonds be- 
yond the critical length correspond to wave vectors at the 
border of the Brillouin zone and ideally should be con- 
nected to the presence of resonances towards the zone 
boundary. Our description however appears to be more 
general as it expresses the same phenomenon in terms of 
the energy release rate and phonon band structure. The 
existence of velocity gaps is shown on the basis of en- 
ergetic arguments. This has the advantage of not being 
built up on the particular rule for breaking bonds, but 
relates the existence of velocity gaps and constant speed 
advance to a more general description of the properties 
of the material. Furthermore, velocity gaps can be read 
directly from the dependence of the strain energy release 
rate on the crack speed, a relationship which we might 
eventually hope to deduce or calculate for real materials. 
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APPENDIX A: DISPERSION RELATIONS FOR 
BULK WAVES 



We need to find the eigenvalues for eq. (2.S). By in- 
spection u k = V A V* is an eigenvector with eigenvalue 



A k = ^[V-V*+V*-V]. 



(Al) 



The remaining two eigenvectors are given by a linear 
combination of V and V*. We then obtain: 

w 2 (W + cV*) = A k (6V + cV*) + 



(a + fj.) si' 



[W* + V*V] • (6V + cV*). (A2) 



We can now multiply eq. ( |A2| ) by V to the left and by 
V* to the right, so as to obtain the following coupled 
equations: 

(e-l)[6rj + c] = {3[2b v + c{l + \ v \ 2 )} 



(e- l)[b + crj* 



p[b(l + M 2 ) + 2ct?*] 



with: e = u; 2 /A k , /3 = (A + fi) SI' V ■ V*/(mA k ) and 
r\ = (V • V)/(V • V*). Solving for eigenvalues leads to: 

e= 1+/3(1± \rj\) with 

and substituting back into the equation we obtain: 

(A + fi) SI' 



to 



A, 



(V ■ V*) ± |V • V| 



(A3) 



which together with flAip completes the dispersion rela- 
tions. 

Having obtained these eigenvalues, it is not difficult to 
find that the corresponding eigenvectors have the form: 



V • V 



(A4) 



APPENDIX B: ADDING DISSIPATION 

There are at least two pratical reasons to add dissi- 
pation in the dynamics. First, the starting point of any 
simulation should be a sample in equilibrium: this can be 
easily obtained by relaxing the lattice to a configuration 
of minimum energy. The second reason is the possibility 
to damp out waves that would otherwise propagate and 
reflect back from the sample borders. 



FIG. 12: Damping rate for different values of 7. Note that 
for k > k* all modes are damped at the same rate q — 7. 



The physical way to introduce dissipation in the model 
for linear elasticity is the introduction of a viscosity term 
in (2.7). However this gives wide dispersion in the damp- 
ing rates of different modes making it inefficient at relax- 
ing solutions or reducing boundary reflections. 

There is a less physical way to obtain dissipation that 
has the advantage of being less wave vector dependent 
and simpler to implement. Starting from eq. ( |2.9| ) a dis- 
sipative term is added so that: 



mii + 2to7U + G(k)u = 



(Bl) 



where: 



G(k)u = Sl'l - fi V • V* + V* • V u 



(A + jti) VV* +V*V 



u 



This approach is clearly unphysical, but useful to reach 
the configuration of minimal energy or to damp out un- 
wanted reflected waves. 

Solutions of the form u(t) = e _9 *u k can be found by 
substitution in (Bl), giving q 2 — 2<yq + G(k)/m = 0, and 
leading to 



Re{q} = Re^7± 



\ 



GOO 



G(k) 



> 



27-^ 
GOO 



G(k) 



2- 



fast 
modes 

slow 
modes. 



GW < , 2 

m. 1 
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In fig. 12 the damping rate Re{<7} is plotted as a func- 
tion of |k|. Curves depend on the value of the dissipa- 
tive parameter 7: in particular, the value A;* (7) sepa- 
rates two regimes. For k > k* all the modes are damped 
with a damping constant q — 7, whereas for k < k* 
fast and slow modes are damped with different magni- 
tudes. The optimal value for 7 is then set by maximising 
the damping rate of the slowest mode, leading us to set 
k* = fc mm = j-, where L is the size of the lattice in fee 
unit cells. The choice is also interesting in that it sets all 
modes to have the same damping rate 



l = c t -. 



(B2) 



The time needed to reduce the potential energy of a factor 
e towards its minimum is then given by: 



At 



h 



2-KC t 



hu 



(B3) 



During the dynamics, reflected waves travel a distance 
of at least L fee unit cells in a time interval of At = L/c t , 
hence when reflected waves reach the crack tip, using the 
optimal value for 7 their intensity is reduced at least by 
a factor e 2jAt = e 2jr ~ 500 for transverse waves, and 

27r £i 

e c i ~ 25 for longitudinal waves, using values as in fig- 
ure H Note however that for the results presented in this 
paper we turned all damping off after initial equilibra- 
tion. 



APPENDIX C: DISPERSION RELATIONS FOR 
SURFACE WAVES 

Consider a solid occupying the portion of space defined 
by y < with the surface y = as its only boundary. 
In the continuum case, surface modes are solutions of 
w 2 Uk + V- er= constrained by the boundary condition 
<7 -n = where n is the normal to the surface. They 
consist of linear combinations of bulk waves with different 
complex k y (that accounts for their damping towards the 
interior of the sample) obeying the bulk waves dispersion 
relations and common u>, k x and k z . 

We can apply the same principles to our discrete model 
considering that the layer of sites corresponding to the 



y = plane has no force acting on it from above. It 
follows therefore: 



Vp CT 2 +Vj- £71= 



(CI) 



where: 



vr = 



?v exp 1 ' 



^2 r v exp 4 ' 

{v\r v>y <0} 



Hence: 



1 Vi = v+ + vr. 

V+(k)^V+(k) = -V 2 -(-k). 



Using eq. ( |2.8| ), the boundary condition (CI) becomes: 
B(k)u = A(V+ • Vi + V+ • V 2 )u+ 

A(ViV+ + V 2 V+) • u + /i(V+Vi + V+V 2 ) • u = 0. 

We can solve this equation by expressing u as a linear 
combination of bulk eigenvectors u = aei + &e 2 + ce 3 
with common u>, k x and k z . For simplicity we can reduce 
to the casec3 k z = 0. 

Bulk dispersion relations for k z = have therefore the 
form u> = f a (k Xl k y ). These can be inverted such that 
k y = /~ (cj, k x ). Hence, surface modes are solutions of 
the following equation: 

/ (B(k) ei 



det 



(B(k)e 2 ) 



V 



B(k)e 3 



ky—ff 3 1 (d),fe x ) / 

provided Im{fp(uj, k x ), / 2 _1 (w, k x ), / 3 _1 (w, k x )} > 0. 
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